function output = arrayEventFinderTest(array, amplitudeThreshold, polarity, durationThreshold)
% arrayEventFinderWithDurThresh
% 
% syntax: 
% output = arrayEventFinderWithDurThresh(array, amplitudeThreshold, polarity, durationThreshold)
%
% array -- e.g. 7000 by 256; rows are data points, columns are traces
% 
% amplitudeThreshold -- set in gui
% 
% polarity -- set in gui; 'negative' or 'positive'
% 
% durationThreshold -- set in gui; in milliseconds
% 
% output --
% output.tracesWithEvents
% output.numberOfEvents
% output.eventOnsetInds
%
% Editing:
% gs april/may 2005 -- modified for mapAnalysis
% gs july 2005 -- re-written from arrayEventFinder
% --------------------------------------------------------------------------------------

if nargin == 0
    data = [0 0 1 1 1 0 1 0 0]'; array = [data data*0 data]; array(8,1) = 1;
    amplitudeThreshold = 0.5; polarity = 'positive'; durationThreshold = 1;
end

[R,C] = size(array);

% take the polarity into account
switch polarity
    case 'positive'
        % no change
    case 'negative'
        % simply flip the polarity to make responses up-going
        array  = -array; 
    otherwise
        disp('Polarity argument must be ''positive'' or ''negative''');
end

for n = 1:C
    data = array(:,n);
    
    % find all the suprathreshold values; e.g. [0 1 1 1 0 1 0 0]
    dataLogical = data > amplitudeThreshold;
    
    % get  indices to all subthreshold regions; e.g. [1 5 7 8]
    dataZero = find(dataLogical == 0);
    
    % differentiate; e.g. [4 2 1]
    diffArrayZero = diff(dataZero);
    
    % subtract 1; e.g. [3 1 0]
    diffArrayZero = diffArrayZero - 1;
    
    % indices of all events; e.g. [1 2]
    indsOfAllEvents = find(diffArrayZero > 0);
    
    % indices of events longer than event duration; e.g. [1]
    indsOfLongEvents = find(diffArrayZero > (durationThreshold));
    allLongEventInds{n} = indsOfLongEvents;
    
    % indices of events shorter than event duration; e.g. [2]
    indsOfShortEvents = setdiff(indsOfAllEvents, indsOfLongEvents);
    
    % number of events
    numberOfEvents(n) = 0 + length(indsOfLongEvents);
    
    % onsets of events
    if indsOfLongEvents
        firstEventInd = indsOfLongEvents(1);
        eventIndsInDataVec = dataZero(firstEventInd + 1) - 1;
        eventDur = diffArrayZero(firstEventInd);
        eventOnsetInds(n) = eventIndsInDataVec + 1 - eventDur;
    else
        eventOnsetInds(n) = NaN;
    end
end

% output:
output.tracesWithEvents = (numberOfEvents > 0)';
output.numberOfEvents = numberOfEvents';
output.eventOnsetInds = eventOnsetInds'; 


% % eliminate short events from traces
% numberOfShortEvents = length(indsOfShortEvents);
% durationsOfShortEvents = diffArrayZero(indsOfShortEvents);
% indsInDataVec = arrayZero(indsOfShortEvents + 1) - 1;
% arrayLogicalProcessed = arrayLogical;
% if ~isempty(indsOfShortEvents)
%     for n = 1:numberOfShortEvents
%         %         indsInDataVec(n) 
%         %         indsInDataVec(n) + durationsOfShortEvents(n)
%         arrayLogicalProcessed(indsInDataVec(n) - durationsOfShortEvents(n) : indsInDataVec(n) ) = 0;
%     end
% end
% [arrayLogical arrayLogicalProcessed];

% % % -------------------- OLD WAY
% % % [R,C] = size(array);
% % % [r,c] = size(amplitudeThreshold);
% % % 
% % % if r*c > 1
% % %     if r > c
% % %         amplitudeThreshold = amplitudeThreshold';
% % %     end
% % %     amplitudeThreshold = repmat(amplitudeThreshold, R, 1);
% % %     array = array-amplitudeThreshold;
% % %     amplitudeThreshold = 0;
% % % end
% % % 
% % % switch polarity
% % % case 'positive'
% % %     [r,c] = find(array >= amplitudeThreshold);
% % % case 'negative'
% % %     [r,c] = find(array <= amplitudeThreshold);
% % % end
% % % 
% % % d = [c r];
% % % dflip = flipdim(d,1);
% % % 
% % % if isempty(d)
% % %     output = [];
% % %     return
% % % end
% % % 
% % % [B,I,J] = unique(dflip(:,1)); % B is a list of all the traces with events
% % % 
% % % output(:,1) = B; % first column: list of all the traces with events
% % % temp = dflip(:,2); 
% % % output(:,2) = temp(I); % second column: time indices of events
% % % output(:,3) = (output(:,1)-1)*size(array,1)+output(:,2);
% % % output(:,4) = array(output(:,3));
% % % events = zeros(C,1);
% % % events(output(:,1)) = output(:,2);
